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ABSTRACT 

We investigate statistics of the decay process in the equal-mass three-body 
problem with randomized initial conditions. Contrary to earlier expectations 
of similarity with "radioactive decay" , the lifetime distributions obtained in 
our numerical experiments turn out to be heavy-tailed, i.e. the tails are not 
exponential, but algebraic. The computed power-law index for the differential 
distribution is within the narrow range, approximately from —1.7 to —1.4, de- 
pending on the virial coefficient. Possible applications of our results to studies 
of the dynamics of triple stars known to be at the edge of disruption are 
considered. 
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1 INTRODUCTION 

Disruption of a three-body gravitational system is an enigmatic dyn amical pr o cess, statis- 
tics of which is mostly unexplored, at long timescales especially. IVaftonenl ( 119881 ) sup- 
posed that the lifetime distribution for a three-body system is an exponentially decay- 
ing function, in anal o gy wi th "radioactive decay". Recent statistical numerical studies by 
Mikkola &: Tanikawal ( 120071 ) of this process in the equal-mass problem revealed new impor- 
tant data on the statistics of the d i srupti on times and raised new questions on the nature of 
this process. Mikkola &: TanikawJ ( 2007 ) have explored a statistics of the disruption times 
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7d in the equal-mass three-body problem. (The disruption time is the system lifetime 
as a bound system.) They computed the disruption times for equal- mass systems with ran- 
domized initial conditions, and fitted the found disruption time distributions by exponential 
functions. 

However, exponenti al decay, in th e asym ptotically long time scale, is in conflict with 



the theoretical result of 



Agekian et al.l ( 119831 ) that the average lifetime of a general i s olated 



three-body system is infinite. What is more, recently it was shown by I Shevchenko! (120091 ) 
on the basis of the Kepler map theory, applied to a hierarchical three-body problem, that 
at the edge of disruption the orbital periods and the size of the orbit of the escaping body 
behave as Levy flights. Due to them, the time decay of the survival probability (the integral 
distribution of lifetimes) is heavy-tailed with the power-law index equal t o —2/3. Combining 



the Kepler map theory and earlier theoretical findings of iHutl ( 119931 ) , IShevchenkol ( 120091 ) 



made a conclusion that the T^ 2 ^ 3 law is expected to be quite universal. 

Here we explore the problem of the lifetime statistics in the three-body problem by 
means of numerical simulations. We pay particular attention to analysis of the tails of the 
lifetime distributions, and find the algebraic behaviour. We compare our results with previous 
numerical work and discuss why the algebraic tails had not been identified earlier. Possible 
applications of our results to studies of the dynamics of triple stars known to be at the edge 
of disruption are considered. 



2 SETTING OF THE PROBLEM AND NUMERICAL EXPERIMENTS 



To investigate disruption of triple systems we use numerical integration of equations of 
motion in the gravitational three-body problem. The equations are as follows: 



d 2 r k 
dt 2 



|^ (* = 1,2,3), 



(1) 



where mk are the masses of the bodies, are the coordinate vectors in the barycentric 
orthogonal coordinate system, t is time, and U is the potential: 



jj _ G f i m 1 m 2 m x m z m 2 m 3 \ 
V 7-ia r X3 r 23 / 



(2) 



where the quantities are the distances between bodies i and j, G is the gravitational 
constant. 

We use the chain-regularization algorithm of iMikkola fc Aarsethl (119931 ) to treat be- 



haviour of the system in the vicinity of close double and triple approaches accurately. Nu- 
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merical integration is performed by the Dormand-Prince integrator (IHairer et al.lll987| ). 
which realizes an explicit Runge-Kutta method of 8-th order with step size control. 

We restrict our present study to the equal-mass three-body problem. The masses of all 
three bodies m^, k = 1, 2, 3 are set to 1. We also choose G = 1. Total energy of the triple 
system E = — 1. 

To s et up the initial cond i tions for the integration, we use an approach similar to that 



used in 



Mikkola fc Tanikawal ( 120071 ) . Namely, we produce random initial data by selecting 



the initial coordinate and velocity components of the bodies from a uniform distribution 
inside a 9-dimensional cube with size of the edge equal to 2. Then we transform these initial 
data to the barycentric system. 

Using scale coefficients, the coordinates and velocities of bodies are transformed to satisfy 
fixed the virial ratio k = T/U, where T and U are the kinetic energy and the potential, 
respectively, of the triple system. The total energy is E = T — U = — 1. The following values 
of the initial virial ratio have been chosen: k = 0, 0.1, 0.3, 0.5, 0.7, 0.9. For each value of the 
virial ratio we study a set of 10 5 triple systems. 

We follow the dynamical evolution of each triple system during the time interval of 10 5 
time uni ts, or until the escape crite rion is satisfied. The escape criterion is the same as that 
used in ( iMikkola fc Tanikawall2007l ): namely, we stop our integration, when a single body is 
moving away on a hyperbolic relative orbit at a distance d > 50 times the current semi-major 
axis of the final binary. The condition for hyperbolicity is 

Eoat = —2~ + ~ G ^^ >0 ' (3) 

where E out is the total energy of the outer binary, M in is the mass of the inner binary, m out 

is the mass of the outer component, V cm is the velocity of centre-of-mass of the inner binary, 
V ont is the velocity of the outer component, d is the distance between centre-of-mass of the 
inner binary and the outer component. 

Using orbital elements we compute the formal time of pericentre passage for this hyper- 
bolic orbit. This time Td is considered to be the moment of the triple system disruption. 
The disruption time is counted from the beginning of computation up to this moment. 



3 RESULTS OF THE NUMERICAL EXPERIMENTS 

The results of integration are qualitatively similar for all values of k in our set, that is why 
the graphs with distributions are presented here only for a single case of k, namely, for k = 0. 
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In Fig. [U we show the integral distribution T a (Td) of the disruption times in decimal 
logarithmic scales. The quantity F a (Td) is the fraction of the lifetimes greater than T d . An 
initial exponential drop (a "bump" in logarithmic scales) and a subsequent power-law decay 
(a straight-line dependence in logarithmic scales) are prominent. In the very tail a smooth 
drop is evident, which is due to the presence of the upper limit on the time of integration 
(equal to 10 5 time units). Owing to this final drop, the F a (T d ) dependence cannot be used 
straightforwardly (by linear fitting of the tail in the logarithmic scales) for reliable estimation 
of the index of the power-law decay. 

We separate the initial exponential drop from the subsequent power-law decay by choos- 
ing the value of T d , at which the distribution becomes algebraic. From Fig. 1 we see that 
this transition value lies between 10 3 and 10 4 . We choose the latter value, so that to be 
completely sure that any contribution of the initial "bump" is minimal. 

Thus the transition value of T d is found from the analysis of the bimodal structure of 
the integral distribution of T d in logarithmic scales. In the subsequent treatment the scales 
are no longer logarithmic. In Fig. [2} the differential distribution of the disruption times is 
shown for the same k = 0. Only the tail (Td > 10 4 ) is presented. The quantity /(Td) is the 
fraction of the lifetimes in a pocket (T d , T d + A), where A = 10 3 . As soon as the distribution 
is differential, the upper limit on the time of integration does not inflict the form of the 
distribution, and so the distribution can be used for fitting. The solid line, drawn in the 



Figure, represents the algebraic fitting. Details of the fitting are given in Table 1. 



kinds of fittings in this paper, we use a non-linear least-squares method ( ILevenberg 



Marquardtl 



' or all 
19441 : 



19631 ) to minimize x 2 > thereby finding the best-fit parameter values and their 
standard errors. 

The fitting here is two-parametric: we use the fitting function 

/(T d ) = AT^, (4) 

where A and ft are free parameters. Only the value of (3 is important for us and is therefore 
reported in Table 1. The correlation coefficient R 2 is also reported. Note that it is very close 

—5 /3 

to 1, i.e. the fitting is very good. As one can see, it turns out to be close to the T d law. 

In Fig. [31 the same data are presented as in Fig. [2j but the distribution is integral. The 
quantity F(T&) is the fraction of the lifetimes less than T d . (Note that the quantity T a (Ti) 
in Fig. 1 is the fraction of the lifetimes greater than Td). The solid line, drawn in the Figure, 
represents the algebraic fitting. A two-parametric fitting cannot be employed in the integral 
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Table 1. The power law index (3 for the tail of the differential distribution of disruption times 



k 





0.1 


0.3 


0.5 


0.7 


0.9 


P 


1.695 
±0.041 


1.639 
±0.032 


1.522 
±0.031 


1.448 
±0.022 


1.586 
±0.028 


1.684 
±0.057 


R 2 


0.967 


0.977 


0.972 


0.984 


0.981 


0.935 



Table 2. The power law index a for the tail of the integral distribution of disruption times 



k 





0.1 


0.3 


0.5 


0.7 


0.9 


a 


0.7073 
±0.0016 


0.6704 
±0.0013 


0.4244 
±0.0022 


0.4555 
±0.0018 


0.6783 
±0.0014 


0.6332 
±0.0041 


R 2 


0.9998 


0.9998 


0.9996 


0.9997 


0.9998 


0.9985 



case, because we do not know beforehand the fraction of the disruption times longer than 
the limiting time of integration (for the differential distribution, this fraction is absorbed in 
the coefficient A). So, the fitting here is three-parametric: we use the fitting function 

F(T d )=A-BTt a , (5) 

where A, B and a are free parameters. Only the value of a is important for us and is 
therefore reported in Table 2. The fitting turns out to be close to T d 2//3 law. This law is 
expected from our results on the differential distribution, because the equality a — (3 — 1 
should hold. In practice this equality is not exact, mostly due to the difference in fitting 
schemes. 

The correlation coefficient in the case of the integral distribution is much closer to 1, 
than in the case of the differential distribution. This means that the estimates of a are 
more reliable than those of (5. This is natural, because the fitting results for the differential 
distribution are sensitive to the choice of the pocket width. 

The results of implementing the fitting procedures at all values of k are collected in 
Tables 1 and 2. As it can be seen from the data in Tables 1 and 2, the tails of distributions 
in all cases are algebraic. This is our first main result, following from the fact that the 
correlation coefficient R 2 in all cases is very close to 1, i.e. the fitting is very good. 

Our second main result is that the computed power-law index for the differential distri- 
butions is in the narrow range, approximately from —1.7 to —1.4, depending on the virial 
coefficient k. The consequences are discussed below. 
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4 DISCUSSION OF THE DISTRIBUTIONS 



What is the nature of the o bserved algebraic decay? One may recall that a similar alge 



braic decay was observed by IShevchenko fc Scholll (119961 . 119971 ) in numerical experiments 



on asteroid dynamics. They showed that the tail of the integral distribution of the time 
intervals T s between jumps of the orbital eccentricity of asteroids in the 3/1 mean mo- 
tion resonance with Jupiter is not exponential, but algebraic: F a oc T~ a with a ~ 1.5—1.7. 
This decay was due to sticking of orbits to chaos borders in the phase space of motion. 
The stickiness effect determines the character of the distribution of Poincare recurrences 



on large timescales: it is algebraic (IChirikov fc 



pioneering work by IChirikov fc Shepelyanskyl ( ]198ll ) , the algebraic decay in the recurrence 



Shepelvanskv 



1981 



1984 ). Starting with the 



st at istics in Hamiltonian syst ems with divided ph a se sp a ce was considered, in par t icular , 



in ( 



Chirikov fc Shepelvanskv 



1984 



Chirikov 



1990 



1999 



Cristadoro fc Ketzmerick 



20081 ). 



Chirikov! ( 119901 ) . using his resonant theory of critical phenomena in Hamiltonian dynamics, 



justified a value of 3/2 for the critical exponent a. 

However, the values of a observed in our numerical experiments (~ 0.4-0.7) differ rad- 
ically from those expected for the stickiness phenomenon (~ 1.5). Therefore, though the 
decay is algebraic, it is not related to the stickiness effect. Then, what is the o rigin of the 
algebraic decay here? On the basis of the Kepler map theory, IShevchenko! (120091 ) considered 
statistics of the disruption and Lyapunov times in a hierarchical three-body problem. It was 
shown that at the edge of disruption the orbital periods and the size of the orbit of the 
escaping body exhibit Levy flights. Due to them, the time decay of the survival probability 
is heavy-tailed with the slope power-law index a equal to 2/3 (while the relation between 
the Lyapunov and disruption tim e s is q u asilinear) . Combini ng the Kepler map theory and 



earlier theoretical findings of iHutl (11993 ). Shevchenko ( 2009 ) concluded that the T d 2 ^ 3 law 



was expected to be quite universal. 



However, neither the T d 2 ^ 3 law, nor even an y other algebraic law, were reported in the 
numerical studies by lMikkola fc Tanikawal (120071 ). We believe t hat the reasons are as f ollows : 



Mikkola & Tanikawa (2007) 



Shevchenko 



solely the initial part of the distribution was considered by 
Beside s, the algebraic fitting functions were not used anyway. As noted by 
J2OO9 ) . the tail of the disruption time distribution in the given problem should be considered 
separately from the initial part, because it corresponds to a different dynamical situation: 
in the beginning, the regime of decay might be exponential; see dynamical analogues in 
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(Shevchenko & Sc 


loll 




L996 


decay in f 


Chirikov 


199< 


)). 



1997 



Shevchenko 



19991 ) and the theory for the exponential 



As follows from Tables 1 and 2, for some values of k, namely intermediate ones k = 0.3 
and 0.5, the power-law index (3 of the differential distribution is equal to 1.45-1.52, while the 
power-law index a of the integral distribution is equal to 0.42-0.46. Thus the values of f3 and 
a are close to 3/2 and 1/2, respectively. This might testify for the presence of the inverse 
square root law F a (T^) oc T d 1 ^ 2 , which is inherent to free diffusi on in the central part of a 
chaot ic layer until th e finite width of th e layer becomes important (I Chirikov fc Shepelyansky 
198ll . p. 9), see also Jshevchenkolll999[ ). 

Concluding, we see that the Kepler map statistics is valid at k ~ and k ~ 1, while, 
hypothetically, the free diffusion seems to dominate at k ~ 0.5. As for the case k ~ 0.5, we 
believe that the free diffusion predominance is merely an effect of insufficient times of inte- 
gration, and expect that an extension of the integration time to greater limits should make 
apparent the asymptotic behaviour typical for the Levy flights in the given problem, when 
appropriately long timescales are achieved; in other words, we expect that the asymptotic 
behaviour with a ~ 2/3 is universal for all values of k. 



5 MULTIPLE STARS: DYNAMICS AT THE EDGE OF STABILITY 



In application to actual multiple stars, our results on the statistics of chaotic decay are 
appropriate to be used in studies of the dynamics of multiple stars known nowadays to be 
near the edge of stability. 

A list of such multiple stars comprises: HP 40887, HP 76644 (A DS 7114= i Uma), 



HP 217675 (o And), HP 222326 (APS 16904) (jOrlov fc Zhuchkov 



20051 ). Probably all these 



systems are quadruple. However, they contain close binary stars, therefore in computer 
simulations we can consider them as triples. The mass ratios in these systems are: 1) 0.65 : 
0.52 : (0.69+0.2); 2) 0.41 : 0.42 : (1.94 + 0.82); 3) 4.2 : 3.4 : (6.8 + 2.3); 4) 1.2 : 1.3 : (0.7+1.0) 
solar masses, respectively. The ratios of the orbital periods of the outer and inner binaries 
are approximate ly 6, 21, 20, 14, respecti vely. If considered as triples, these systems have a 
weak hierarchy (jOrlov fc Zhuchkov! 120051 ): the distant component is not very much far from 
the inner pair (with respect to the size of the latter), and its perturbation on the motion of 
the latter is not negligible. The mutual perturbations can be strong enough to modify the 
hierarchy and, finally, lead to disruption of the triple. Thus the property of weak hierarchy 
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implies that these systems are near t he edge of the stability region in the phase space of the 
motion; see (jOrlov fc Zhuchkovll2005l ) for details. Though different from the equal-mass case, 
the dynamical behaviour of these systems should be as well dominated by Levy flights, if 
they are close enough to disruption. This proximity should be further analyzed and verified. 

Of course, a direct observation of Levy flights in the dynamics of such systems would 
require unrealistic long times of observation. However, when the current coordinates and 
velocities of the components are determined observationally with sufficient accuracy, the 
Levy flights can be studied and analyzed in numerical simulations of the future dynamical 
history of the triples, and estimates of their lifetimes can be obtained. 



6 CONCLUSIONS 



We have investigated statistics of the disruption times Ta in the equal-mass three-body 
problem. The statistics has been explored in massive numerical simulations with randomized 
initial conditions. The distributions of the disruption times have turned out to be heavy- 
tailed with the power-law index being within narrow range near the value of —5/3 (for the 
differential distribution). The observed range is from —1.7 to —1.4; the value of the slope 
index slightly depends on the virial coefficient. 

Our result is in conflict with earlier expectations and analogies with "ra dioactive decay" , 



Mikkola fc Tanikawa 



which implied exponential tails, and also in conflict with numerical results ofy 
(120071 ). where ex ponential decay was found. W e explain the divergence with the latter results 
by the fact that iMikkola fc Tanikawa! (120071 ) fitted the distribution as a whole (and what 
is more, taken at insufficiently long timescale), while the behaviour of the tails should be 
treated separately from the initial exponential drop. 

On t he other hand, pur fin ding of the algebraic decay is in agreement with the theoretical 



result of 



Agekian et al. 



( 119831 ). stating that the average lifetime of a general is olated three- 



body system is infinite. Moreover, it is in agreement with the Kepler map theory (jShevchenko 
20091 ). predicting just the same value of the slope index (~ —5/3) that we have observed in 



the majority of our numerical experiments. This makes evident the existence of Levy flights 
in the process of decay. 

At intermediate values of the virial coefficient, k ~ 0.5, we have observed the slope index 
~ —1.5 (for the differential distribution). This is typical for free diffusion in the central 
part of a chaotic layer until the finite width of the layer becomes important. We believe 
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that the free diffusion predominance is merely an effect of insufficient times of integration. 
If one extends the integration time to sufficient limits, the asymptotic behaviour with the 
slope index ~ —5/3, characteristic of the Levy flights in the given problem, should become 
apparent. 

Several actual multiple stars are known nowadays to be near the edge of the stability 
region in the phase space of the motion, in particular, HP 40887, HP 76644 ( ADS 7114= 
l Uma), HP 217675 (o And), HP 222326 (APS 16904) Jorlov & Zhuchkov! 120051 ). Numerical 
integrations, based on initial data from further high-precision astrometric observations of 
these systems, would allow one to reveal long-term evolution features, in particular, the 
theoretically predicted Levy flights in the orbital periods and the geometrical sizes of the 
systems. 
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Figure 1. The integral distribution of the disruption times in logarithmic scales, k = 0. 
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Figure 2. The differential distribution of the disruption times for k = 0. The solid line represents the algebraic fitting, which 

— 3/2 

turns out to be close to T d law. 
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Figure 3. The same as in Fig. [2] but the distribution is integral. The solid line represents the algebraic fitting 



